Hybrid elastic/discrete-particle approach to biomembrane dynamics with application to the 

mobility of curved integral membrane proteins 



Ali Naji, 1 Paul J. Atzberger, 2 and Frank L. H. Brown 1 

'Department of Chemistry and Biochemistry & Department of Physics 
2 Department of Mathematics, University of California, Santa Barbara, CA 93106, USA 

We introduce a simulation strategy to consistently couple continuum biomembrane dynamics to the motion 
of discrete biological macromolecules residing within or on the membrane. The methodology is used to study 
the diffusion of integral membrane proteins that impart a curvature on the bilayer surrounding them. Such 
proteins exhibit a substantial reduction in diffusion coefficient relative to "flat" proteins; this effect is explained 
by elementary hydrodynamic considerations. 
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Lipid-bilayer membranes are among the most important 
and most versatile components of biological cells HI Q]; 
biomembranes protect cells from their surroundings, provide 
a means to compartmentalize subcellular structures (and the 
functions of these structures) and act as a scaffolding for 
countless biochemical reactions involving membrane associ- 
ated proteins. Our conceptual picture of biological mem- 
branes as a "two-dimensional oriented solution of integral 
proteins ... in the viscous phospholipid bilayer" iQl was popu- 
larized well over thirty years ago. Quantitative physical mod- 
els for the energetics J3l and dynamics ^ associated with 
shape fluctuations of homogeneous fluid membranes and for 
the lateral diffusion coefficient of integral membrane proteins 
within a flat bilayer were developed shortly thereafter and 
still find widespread use up to this day. 

Interestingly, the coupling of protein diffusion to the shape 
of the membrane surface has become a subject of study only 
relatively recently (see 0, H Ell and references within). One 
well studied consequence of membrane shape fluctuations is 
that a protein must travel a longer distance between two points 
in 3D space if the paths connecting these points are con- 
strained to lie on a rough surface as opposed to a flat plane. 
This purely geometric effect is expected to have practical ex- 
perimental implications; measurements that capture a 2D pro- 
jection of the true motion over the membrane surface will 
infer diffusion coefficients of diminished magnitude relative 
to the intrinsic lateral diffusion locally tangent to the bilayer 
surface J3, H 0, OH HI- Beyond this generic effect, which 



should apply to anything moving on the membrane surface in 
relatively passive fashion (lipids, proteins, choleseterol, etc.), 
certain membrane associated proteins effect shape changes in 
the bilayer. Specific examples include the SERCla calcium 
pum p Il2ll and BAR (Bin, Amphiphysin, Rvs) domain dimers 
1 13]. Direct structural evidence from X-ray crystallography 



112111311, experimental studies of vesicle topology at the mi- 
cron scale flii [l4ll . and atomically detailed simulations 1 151 



all indicate the ability of these proteins to drive membrane 
curvature (see Fig. [T). 

The diffusion of membrane proteins with intrinsic curva- 
ture is more complex than the diffusion of a relatively pas- 
sive spectator and remains incompletely explored in the liter- 



ature. Although stochastic differential equations coupling the 
lateral motion of curved proteins to thermal shape fluctuations 
of a continuous elastic bilayer have been proposed |n, Trill , 
these equations have only been analyzed under the simplify- 
ing assumption that the protein does not affect the shape of the 
membrane surface. Under such an approximation, it was pre- 
dicted both analytically ll 111 and numerically tfl6f | that curved 
proteins are expected to diffuse more rapidly than flat ones. 
The full numerical analysis provided in the present work sug- 
gests exactly the opposite effect-protein curvature decreases 
lateral mobility across the bilayer. The disagreement with pre- 
vious work is attributable to the fact that the protein's influ- 
ence on bilayer shape is of primary importance (see Fig. Q]) 
and cannot be ignored. 

Our starting point is the Monge-gauge Helfrich Hamilto- 
nian J3l for energetics of a homogeneous membrane surface 
under conditions of vanishing tension 



H = \ [ dx [A m {V 2 h) 2 + 2K' m 0(h)] . 
1 Ja ± 



(1) 



where /i(x) = h(x, y) describes the local membrane displace- 
ment from a flat reference plane at z = (see Fig. [T). Here, 
Q{h) = {d xx hd yy h — d xy hd xy h) is the Gaussian curvature 
and K rn and K' m are the membrane bending modulus and 
saddle-splay modulus, respectively. The integration region, 
A± = L 2 , is always taken to be a square box with peri- 
odic boundary conditions assumed. We treat membrane pro- 
teins as localized regions of enhanced rigidity within the bi- 
layer. A single protein centered within the bilayer at position 
(r(i), h(r(t))) = (x(t), y(t), h(x(t), y(t))) is thus assumed 
to modify the Hamiltonian as H = Hq + Hint with 

H int = I [ dxG p (x-r(t)) [K p (V 2 /i - 2C P ) 2 (2) 
1 Ja ± 

-K m (W 2 h) 2 + 2(K' p -K' m )g(h)]. 

K p and K' p are the protein bending and saddle-splay moduli 
and C p is the spontaneous curvature associated with protein 
shape. The protein shape function G p (x — r(i)) describes the 
envelope of protein influence over bilayer elastic properties; 
the specific function chosen will be discussed in detail below. 
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FIG. 1: (Color online) An appropriately shaped protein will tend 
to distort the local shape of the bilayer into a bent configuration. In 
our simulations, the protein's effective size is defined via the enve- 
lope function in Eq. l|5} (displayed next to the cartoon, top). The 
distortion is most readily seen by minimizing TL for the composite 
protein-bilayer system (middle), however the protein's influence is 
strong enough to maintain visibly apparent perturbations, even in the 
presence of fluctuations (bottom). See Table|I]for parameters. 

Using a Fourier representation /i(x) = p^ q 'iqe" , ' x , 
membrane dynamics may be cast as a set of coupled Langevin 
equations for the individual Fourier modes lfl7i 18ll 



M*) = A q F q (i) + x /fc B TL 2 A q £ q (i), 



(3) 



where Fq(t) is the Fourier-transform of the force per unit area 
F(x,t) = —8TL/8h{x.,t), A q = 1/(4?7<j) corresponds to the 
Oseen hydrodynamic kernel A(x) = l/(87r7y|x|), and £ q (i) is 
a Gaussian white noise with (£ q (<)) = and (£ q (£) £ q ' (t 1 )) = 
2<5 qj _ q ' 5(t — t'), to ensure satisfaction of the fluctuation- 
dissipation theorem. 

The protein's position may similarly be described in terms 
of a Langevin equation that implicitly enforces protein local- 
ization to the membrane surface JH ■ The two independent 
variables describing this motion are the components of r(t) 
(with i, j = 1,2 and summation convention assumed) 



h(t) =D oVi + v^Ty Vj (t) + -^(g-'hf: 



JIJJ3- 



(4) 



Here (g 

sor, g = 1 



Sij — dihdjh/g is the inverse metric ten- 
(V/i) 2 is the determinant of the metric 
gij = 5ij + dihdjh and we have defined di = d/dxi. 
Equation (0]i introduces geometric factors including Vi = 
- [{9~ l )]k djd k h] dih/g and = <5 y - dihdjh/ (g + y/g) 
(i.e., the square-root of the inverse metric tensor (<? _1 )y = 
TikTjk). The Gaussian white-noise rji(t) with (r)i(t)) = 
and (rji(t) rjj(t')) = Sij 8{t — t') guarantees that a tracer 
particle (defined by 7ij nt = 0) undergoes a curvilinear ran- 
dom walk over the membrane surface with diffusion coeffi- 



cient Dq- The last term in Eq. (|4|i reflects the effect of the 
interaction-induced force fa = —dHi n t/dxi. 

Equations (H}-© specify the stochastic thermal evolution 
for a single curved protein coupled to an elastic membrane. 
As noted above, equations very similar to these have been 
proposed previously 1 11, TfJ, but have only been analyzed by 
neglecting the contribution of 7ii nt to F q while maintaining 
its contribution to /j. To avoid this uncontrolled approxima- 
tion, it is necessary to introduce a numerical algorithm that 
can consistently couple protein position r(t) to membrane un- 
dulations /i(x, t). 

For both physical and numerical purposes, we must truncate 
the membrane modes at some short-distance scale a, which 
can for example be taken to represent the typical molecular 
(lipid) size or bilayer thickness. We thus limit the Fourier 
modes appearing in Eq. ([3]l to q = (q x ,q y ) = (27m, 2nm)/L 
where M = L/a with integer n,m in the range — M/2 < 
n,m < M/2. In principle, this reduced set of modes de- 
scribes a fully continuous membrane height profile via /i(x) = 
■jrj hq e lq x at any given point x. However, it is computa- 
tionally advantageous to explicitly track the membrane height 
only over the discrete M x M real-space lattice (defined at 
positions x Q = (p, q)a with integer < p, q < M) conjugate 
to the chosen /i q 's via Fast-Fourier Transformation IU9I1 . 

The interaction between a fully continuous variable de- 
scribing protein position r(i) and a discrete representation 
of the membrane height field {h{x. a ,t}} poses certain chal- 
lenges. A minor issue is that Eq. requires the shape of 
the membrane surface over the entire x, y plane and not just 
at the lattice sites {x Q }. This problem is readily handled via 
linear interpolation to obtain h(x, t) and the required deriva- 
tives at arbitrary x from the corresponding neighboring lat- 
tice values [JsJ] . A more complex problem involves the dy- 
namics of the hq(t)'s in Eq. (|3). The forces in this equa- 
tion include contributions due to the coupling between pro- 
tein and bilayer from Eq. (0. The envelope function G p (x) 
reflects protein size and is quite localized in real space; the 
natural way to deal with Eq. (|2]i (and the related force expres- 
sions) is to approximate the integral by simple quadrature, i.e. 
/ dxG p (x - r(t))^(x) w a 2 £ Q G p (x Q - r(t))JF(x a ) for 
arbitrary function T. To define a numerical scheme that is 
both efficient and accurate, the specific functional form cho- 
sen for Gn is critical. Naive continuous choices like 2D Gaus- 



sians IU61I18I1 lead to an effective normalization (as computed 
by quadrature) that varies with the offset between the enve- 
lope center and the discrete lattice. Piecewise linear forms 
for G p can be defined that suffer no such normalization is- 
sue, but such functions lead to discontinuous derivatives as 
the protein-lattice offset changes. Both scenarios are unac- 
ceptable as these numerical issues lead to a breaking of the 
homogeneity of the membrane surface; the protein will tend 
to favor (or disfavor) lattice sites over other regions of space. 

Our numerical description of coupled membrane-protein 
dynamics shares features with the Immersed Boundary for- 
mulation of hydrodynamics ll20ll . In that work, a series of 
envelope functions are introduced that are continuous, local- 



3 



parameter 


box 


lattice 




protein 


temperature 


bare protein 


bilayer 


protein 


saddle-splay 


protein 


solvent 




dimension 


spacing 




area 




diffusion 


bending 


bending 


moduli 


spontaneous 


viscosity 














coefficient 


modulus 


modulus 




curvature 




symbol 


L 


a 


M = i 

a 




T 


Do 


K m 


K P 


K' 


c P 


V 


value 


250nm 


2.5nm 


100 


lOOnm 2 


300K 


5.0^- 


5k B T 


40k B T 


-Km(p) 


O.lnrrr 1 


»7w=0.001Pa ■ s 



TABLE I: Default parameter values used in the simulations. The protein area is chosen as A p — (4a) 2 (giving a protein diameter of a p ~ 
lOnm), the Do value qualitatively reflects the motion of band 3 protein dimers on the surface of human red blood cells fiUl . and r]^ stands for 
the viscosity of water. For simplicity we assume that both bilayer and protein saddle-splay moduli are equal in magnitude but opposite in sign 
to the corresponding bending moduli. 




ized in space and strictly preserve normalization as evalu- 
ated by quadrature. For our purposes, we take G p (x) = 
{A p /a 2 )cp{x/a)cp{y/a) with Q 

-2 < u < -1, 
-1 < u < 0, 

< u < 1, 

1 < u < 2, 

(5) 

and zero for all other u values (see Fig. [TJ for a plot). A p 
defines an effective protein area and the envelope function is 
non-vanishing over a total of 64 lattice sites. In order to simu- 
late the dynamics of the system, we evolve Eqs. ([3]) and in 
time via the Euler-Maruyuma method 12111 . The resulting al- 
gorithm is essentially an application of "Brownian dynamics 
with hydrodynamic interactions" 1I22I1 applied to membrane 
shape fluctuations and protein motion. Adopting the envelope 
function defined in Eq. © allows a seamless melding of the 
approaches introduced in Refs. JSHHl (a detailed descrip- 
tion of our algorithm will be provided in a future publication). 

The protein's influence on the membrane is most easily seen 
in the absence of thermal shape fluctuations of the membrane, 
but for sufficiently large C p and K p values the effect of the 
protein remains clearly visible despite thermal fluctuations of 
the membrane surface (see Fig. [TJ|. The average distortion of 
the membrane surrounding the protein is sufficient to signif- 
icantly slow protein motion in all cases we have studied (see 
Fig. This slowing derives from two effects. First, the pro- 
tein tends to trap itself in the deformation created by its own 
perturbation to the bilayer. The energy-minimized configura- 
tion displayed in Fig. [TJplaces the protein at the bottom of a 
curved valley; attempted diffusion up the walls of this valley 
is hindered by the interaction-induced force in Eq. (0|. Sec- 
ond, x, y translation of the protein is accompanied by trans- 
lation of the membrane deformation surrounding this protein. 
The hydrodynamic effects included in Eq. <f3j dictate that the 
translation of such a deformation is resisted by the viscous 
drag of the medium surrounding the bilayer. This drag acts 
in addition to the usual quasi-2D drag incorporated within Dq 
and slows protein motion. These effects are most pronounced 
when shape fluctuations of the bilayer are neglected by setting 
T = in Eq. ([3]). Increasing the local deformation around the 
protein or the solvent viscosity both decrease D. In Fig. [2] we 
display three means to control this effect. Increasing C p drives 



large deformations from a flat plane for any finite K p . Larger 
values of K p will tend to increase this effect up until the point 
where the rigidity of the protein becomes effectively infinite 
and the response to the protein saturates. The hydrodynamic 
drag in our model is controlled via rj; increases in -q reduce D 
as effectively as do perturbations to membrane shape. 

Within the asymmetric coupling approximation, which ig- 
nores the influence of the protein on membrane shape, it is 
predicted that bilayer shape fluctuations will enhance curved 
protein mobility UH EH] (i.e. D > D , insets Fig. |2|. Al- 
though this effect can be seen in our simulations, the enhance- 
ment is slight (relative to the similar effect within the afore- 
mentioned approximation scheme) and can not overcome the 
dominant slowing caused by the protein's distortion of the bi- 
layer (compare open circles and triangles, Fig. [2}. We find 
D < Dq for all cases studied, which represents a qualitative 
departure from earlier predictions. 

We may approximately account for the viscous drag ef- 
fect discussed above by invoking an adiabatic approximation 
and assuming that the energy-minimized membrane distortion 
profile (denoted by h(x, t)) instantaneously tracks protein po- 
sition. Hence, 



h(x,t) = -v ■ V/i(x,t), 



(6) 



for a protein moving at constant lateral velocity v. The power 
dissipated by viscous losses in the medium may be calcu- 
lated from P = TjEqM^-qW' where *q(*) follows 
from Eq. (f3]) as F q = /i q /A q . Thus, by using Eq. (O, the 
power loss may be written as P = X) q ( v ' l) 2 l^qlV-^-q — 
|v| 2 //Xdcf with /idcf being the effective mobility of the de- 
formation. The effective protein diffusion coefficient follows 
from D def = HdcikuT and D^ 1 = (D Q 1 + D^ f ) to give 

1 ^ 1 277 
D ~ D k B TL 2 



(7) 



This expression depends only on the minimized deformation 
profile of the membrane at fixed protein position and is read- 
ily calculated (shown as solid squares in Fig. |2j. Although 
the approximation is imperfect due to neglect of membrane 
fluctuations and the self-trapping effect discussed above, the 
adiabatic results serve as a reliable estimator of the observed 
trends for the full simulation. 

We are not aware of experimental studies that specifically 
investigate the role of protein curvature on self-diffusion, but 
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FIG. 2: Projected diffusion coefficients for protein motion on an elas- 
tic membrane. In each panel, one of the three parameters K p (pro- 
tein bending modulus), C p (protein spontaneous curvature), and r\ 
(medium viscosity) is varied as indicated while holding the remain- 
ing physical properties at the default values indicated in TableU Cir- 
cles indicate the results of full simulations including thermal motion 
of both the bilayer and the protein. The triangles indicate results ob- 
tained by turning off all bilayer shape fluctuations (i.e. setting T — 
for membrane undulations modes) and the solid squares indicate the 
results of the adiabatic theory discussed in the text. Insets show sim- 
ulation results obtained within the asymmetric coupling approxima- 
tion. Errorbars are approximately the size of the symbols. Diffusion 
coefficients are calculated as D = [r(t) — r(0)] 2 /(4t). The simula- 
tions run for ~ 5 x 10 7 time steps of size At ~ 0.01ns, ensuring both 
numerical accuracy and statistically reliable results. The protein's 
root-mean-square displacements over the course of the simulations 
are approximately lOOnm (i.e. 20 times the protein's radius). 



do note that recent experiments J23ll show deviations from the 
standard theory |0] used to predict membrane-protein mobil- 
ity. The effect described here may prove to be important in 
describing these deviations for certain proteins. The solvent 
viscosity dependence we find is at odds with the weak (log- 
arithmic) dependence expected for flat proteins |@1 and pro- 
vides a concrete means to verify our predictions experimen- 
tally. 
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